Numerical simulation of structural behaviors using a meshfree-enriched finite element method

ABSTRACT

System, method and software product for numerically simulating structural behaviors of an engineering product in compressible and near-incomprssible region is disclosed. Meshfree enriched finite element method (ME-FEM) is used for such numerical simulation. ME-FEM requires an engineering product be represented by a FEM model comprising a plurality of finite elements. Finite elements used in the ME-FEM are generally low-order finite elements. Each of the finite elements in the FEM model is enriched by at least one meshfree enriched (ME) node located within the element&#39;s domain. Each ME node has additional degrees-of-freedom for the element it belongs independent from those of the corner nodes. A displacement based first-order convex meshfree approximation is applied to the ME node. The convex meshfree approximation has Knonecker-delta property at the element&#39;s boundary. The gradient matrix of ME-FEM element satisfies integration constraint. ME-FEM interpolation is an element-wise meshfree interpolation that is discrete divergence-free at the incompressible limit.

FIELD OF THE INVENTION

The present invention generally relates to computer-aided engineering(CAE) analysis of engineering products or components, more particularlyto systems and methods for numerical simulation of structural behaviorsusing a meshfree-enriched finite element method (ME-FEM).

BACKGROUND OF THE INVENTION

Traditionally, FEM or Finite Element Analysis (FEA) is one of the mostpopular computer aided engineering tool for engineers and scientists tomodel and solve engineering problems relating to complex systems such asthree-dimensional non-linear structural design and analysis. FEA derivesits name from the manner in which the geometry of the object underconsideration is specified. With the advent of the modern digitalcomputer, FEA has been implemented as FEA software. Basically, the FEAsoftware is provided with a model of the geometric description and theassociated material properties at each point within the model. In thismodel, the geometry of the system under analysis is represented bysolids, shells and beams of various sizes, which are called elements.The vertices of the elements are referred to as nodes. The model iscomprised of a finite number of elements, which are assigned a materialname to associate with material properties. The model thus representsthe physical space occupied by the object under analysis along with itsimmediate surroundings. The FEA software then refers to a table in whichthe properties (e.g., stress-strain constitutive equation, Young'smodulus, Poisson's ratio, thermo-conductivity) of each material type aretabulated. Additionally, the conditions at the boundary of the object(i.e., loadings, physical constraints, etc.) are specified. In thisfashion a model of the object and its environment is created.

However, FEM in certain circumstances has its drawbacks. For instance,FEM may not be advantageous due to numerical compatibility condition isnot the same as the physical compatibility condition of a continuum. Inparticular, in Lagrangian type of computations, one may experience meshdistortion, which can either end the computation altogether or result indramatic deterioration of accuracy. In addition, the FEM often requiresa very fine mesh in problems with high gradients or a distinct localcharacter, which can be computationally expensive. In a procedure calledArbitrary Lagrangian Eulerian (ALE) formulations, its objective is tomove the mesh independently of the material so that the mesh distortioncan be minimized. But the distortion in general still creates severalnumerical errors for large strain and/or high speed impact simulations.Furthermore, a mesh may carry inherent bias in numerical simulations,and its presence becomes a nuisance in computations. An example is thestrain localization problem, which is notorious for its mesh alignmentsensitivity. Therefore it is computationally efficacious to discretize acontinuum by a set of nodal points without mesh constraints.

As a result, in recent years, meshfree method has also been used tosolve these drawbacks suffered in FEM. However, meshfree method suffersits own problems. For example, the computational costs are generally toohigh comparing to that of FEM. It is very difficult and computationalexpensive to impose essential boundary conditions in meshfree method.

Neither FEM nor meshfree method can numerically predict every possiblestructural behaviors under various environmental loads, for example, anumerical phenomenon called volumetric locking for simulating structuralbehaviors in compressible and/or near-incomprssible regions (e.g.,rubber, material under large plastic deformation). FIG. 1A depicts asimplified two-dimensional example illustrating the volumetric lockingeffect. Element 102 is fixed in two of its four sides. As a result, node112 is fixed in space thus not moving in numerical simulation. However,this does not represent a physical reality.

Volumetric locking can be described using a set of mathematicalderivation as follows:

Consider a plane-homogeneous isotropic linear material body whichoccupies a bounded polygonal domain Ω in

² with boundary Γ=∂Ω. The problem is restricted to small deformation andthe infinitesimal strain tensor g is defined as a function of thedisplacement u by

$\begin{matrix}{{ɛ(u)} = {\frac{1}{2}\left( {{\nabla u} + {\nabla u^{T}}} \right)}} & (1)\end{matrix}$For a prescribed body force f, the governing equilibrium equation in Ωreads∇·σ+f=0  (2)where σ is the symmetric Cauchy stress tensor. With the fourth-orderelasticity tensor denoted by C, the constitutive equation is given byσ=Cε=2με+λ(trε)I in Ω  (3)where I is the identity tensor of order four. Symbols μ and λ are Laméconstants which are related to the Young's modulus E and Poisson ratio vby

$\begin{matrix}{{\mu = \frac{E}{2\left( {1 + v} \right)}},{\lambda = \frac{vE}{\left( {1 + v} \right)\left( {1 - {2v}} \right)}}} & (4)\end{matrix}$

For simplicity, the displacement is assumed to satisfy the homogeneousDirichlet boundary condition u=0 on Γ_(D) and the traction t isdistributed over the Neumann boundary Γ_(N) with Γ=Γ_(D)∪Γ_(N) andΓ_(D)∪Γ_(N)=0. The variational form of this problem is to find thedisplacement uεV such thatA(u,v)=l(v)∀vεV  (5)where the space V=H₀ ¹(Ω) consists of functions in Sobolev space H¹(Ω)which vanish on the boundary and is defined byV(Ω)={v:vεH ¹ ,v=0 on Γ_(D)}  (6)

The bilinear form A(.,.) and linear functional l(.) in Eq. (5) aredefined by

$\begin{matrix}\begin{matrix}{{A:{VxV}}->} & {{A\left( {u,v} \right)} = {{2\;\mu{\int_{\Omega}^{\;}{ɛ(u)}}}:{{{ɛ(v)}\ {\mathbb{d}\Omega}} + {\lambda{\int_{\Omega}^{\;}{\left( {\nabla{\cdot u}} \right)\left( {\nabla{\cdot v}} \right)\ {\mathbb{d}\Omega}}}}}}}\end{matrix} & (7) \\\begin{matrix}{{l:V}->\bullet} & {{l(v)} = {{\int_{\Omega}^{\;}{{f \cdot v}\;{\mathbb{d}\;\Omega}}} + {\int_{\Gamma_{N}}^{\;}{{t \cdot v}\ {\mathbb{d}\Gamma}}}}}\end{matrix} & (8)\end{matrix}$

The bilinear form A(.,.) in the above equation is symmetric, V-ellipticand continuous on V. By Korn's inequality and Lax-Milgram theorem, thereexists a unique solution to the problem.

The standard Galerkin method is then formulated on a finite dimensionalsubspace V^(h)⊂V employing the variational formulation of Eq. (5) tofind u^(h)εV^(h) such thatA(u ^(h) ,v ^(h))=l(v ^(h))∀v ^(h) εV ^(h)  (9)

The second term on the right hand side (RHS) of Eq. (7) resembles thepenalty term similar to the classical penalty method in the Stokesproblems. For the energy in Eq. (7) to remain finite as λ→∞ (or v→0.5),the following constraint must be enforced∇·u=0 for uεV  (10)Similarly in the finite dimensional space, we have∇·u ^(h)=0 for u ^(h) εV ^(h)  (11)

The solution of Eq. (9) using a low-order approximation for displacementfield in general does not satisfy the constraint in Eq. (11) and failsto converge uniformly as v→0.5. This is the volumetric locking in thenear-incompressible problem. In other words, volumetric locking occurswhen the approximation space V^(h) is not rich enough for theapproximation to verify the divergence-free condition ∇·u^(h)=0. Shownin FIG. 1B are eight non-physical modes for 3×3 discretization usingEFG—a prior art meshfree procedure.

An eigenvalue analysis using the left hand side (LHS) of Eq. (9) can beemployed to further reveal the locking phenomena in the low-orderdisplacement-based finite element and meshfree methods. The eigenvalueanalysis of one rectangular quadrilateral bilinear (Q1) finite elementwith 2×2 Gauss integration points is analyzed first and two non-physicallocking modes are obtained and shown in FIGS. 1C and 1D. The eigenvaluesgrow unbounded as v approaches to 0.5, which is not expected physically.Accordingly, eight non-physical locking modes are obtained from theeigenvalue analysis using four rectangular Q1 elements. Similar to theresults in the four Q1 finite element model, the same number ofnon-physical locking modes are generated by Element-Free Galerkin (EFG)method as displayed in FIG. 1B using the first-order approximation withnormalized nodal support equal to 1.5. The increase of nodal support andorder of Gaussian quadrature rule does not completely eliminate thevolumetric locking in EFG method. The number of non-physical modesremains unchanged.

Based on the aforementioned drawbacks, shortcomings and problems, itwould, therefore, be desirable to have an improved computer-implementedmethod for numerically simulating structural behaviors in compressibleand/or near-incompressible regions.

BRIEF SUMMARY OF THE INVENTION

This section is for the purpose of summarizing some aspects of thepresent invention and to briefly introduce some preferred embodiments.Simplifications or omissions in this section as well as the abstract andthe title herein may be made to avoid obscuring the purpose of thesection. Such simplifications or omissions are not intended to limit thescope of the present invention.

The present invention discloses a system, method and software productfor numerically simulating structural behaviors of an engineeringproduct in compressible and/or near-incomprssible region.

According to one aspect, meshfree enriched finite element method(ME-FEM) is used for such numerical simulation. ME-FEM requires anengineering product be represented by a FEM model comprising a pluralityof finite elements (e.g., plate, solid, etc.). The finite elements usedin the ME-FEM are generally low-order finite elements (i.e., finiteelements having only corner nodes). Each of the finite elements in theFEM model is enriched by at least one added node referred to as meshfreenode or meshfree enriched (ME) node located within the element's domain(i.e., within the outer boundary of the element).

According to another aspect, each ME node has additionaldegrees-of-freedom for the element it belongs independent from those ofthe corner nodes. A particular procedure referred to as displacementbased first-order convex meshfree approximation is applied to the MEnode. The convex meshfree approximation has Knonecker-delta property atthe element's boundary and is constructed using a special meshfreetreatment referred to as generalized meshfree (GMF) procedure. Thegradient matrix of ME-FEM element satisfies integration constraint. TheME-FEM interpolation is an element-wise meshfree interpolation that isdiscrete divergence-free at the incompressible limit. Divergence-freeproperty is provided by meshfree enrichment in quadrilateral elements.For triangular and tetrahedral elements, divergence-free property isprovided by a combination of meshfree enrichment and strain smoothingprocedure.

According to another aspect, the displacement based first-order convexmeshfree approximation uses at basis function derived from anexponential, an inverse tangent, a hyperbolic tangent or Renyi function.

Optionally in another aspect, the finite elements can be divided intofirst and second groups. Only the finite elements belong to the firstgroup are enriched (i.e., adding one or more ME nodes). In other words,ME-FEM can include regular finite elements (i.e., second group) andmeshfree enriched elements (i.e., first group).

According to yet another aspect, the ME-FEM element is configured toavoid possible pressure oscillation in simulations of structuralbehaviors in near-incompressible region. To achieve this goal, anarea-weighted integration is incorporated with the divergence-freeME-FEM interpolation, which provides smoothing effect on strains andpressures. The area-weighted integration further providesdivergence-free property particularly in triangular and tetrahedralelements.

According to still another aspect, a modified Hu-Washizu variationalprinciple is used for creating discrete equations with the smoothedstrain field.

As a result of using ME-FEM elements, numerical simulations ofstructural behaviors in compressible and/or near-incompressible regioncan be performed with reasonable accuracy for assisting user (i.e.,engineers, scientists, etc.) to improve design of an engineering product(i.e., car, airplane, etc.).

Other objects, features, and advantages of the present invention willbecome apparent upon examining the following detailed description of anembodiment thereof, taken in conjunction with the attached drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features, aspects, and advantages of the presentinvention will be better understood with regard to the followingdescription, appended claims, and accompanying drawings as follows:

FIG. 1A is a diagram showing a finite element to demonstrate volumetriclocking;

FIG. 1B is a diagram showing various non-physical modes for 3×3discretization;

FIGS. 1C and 1D are diagrams showing two non-physical modes for Q1bilinear finite element;

FIG. 2A is a diagram showing an exemplary domain is divided ortriangulated into a number of finite elements;

FIG. 2B is a diagram showing an exemplary isoparametric mapping schemeused in a 5-noded meshfree-enriched finite element (quadrilateralelement) in accordance with one embodiment of the present invention;

FIG. 2C is a diagram showing an exemplary isoparametric mapping schemeused in a 4-noded meshfree-enriched finite element (triangular element)in accordance with one embodiment of the present invention;

FIG. 2D is a diagram showing an exemplary isoparametric mapping schemeused in a 5-noded meshfree-enriched finite element (tetrahedral element)in accordance with one embodiment of the present invention;

FIGS. 3A and 3B are diagrams showing shape functions for an exemplarymeshfree-enriched finite element in accordance with one embodiment ofthe present invention;

FIG. 4 is a diagram showing relationship between an interior node withits neighbor elements in accordance with an embodiment of the presentinvention;

FIGS. 5A and 5B show a non-physical mode using the original meshfreemethod and a realistic mode using ME-FEM, according to an embodiment ofthe present invention;

FIGS. 6A and 6B are diagrams showing exemplary area-weighted integrationprocedures used in an embodiment of the present invention;

FIG. 7 is a diagram showing eight locking-free modes in four ME-FEMelements including area-weighted integration;

FIG. 8 is a flowchart illustrating an exemplary process of numericalsimulation of structural behaviors in compressible ornear-incompressible region using ME-FEM in accordance with oneembodiment of the present invention; and

FIG. 9 is a function diagram showing salient components of an exemplarycomputer system, in which an embodiment of the present invention may beimplemented.

DETAILED DESCRIPTION

In the following description, numerous specific details are set forth inorder to provide a thorough understanding of the present invention.However, it will become obvious to those skilled in the art that thepresent invention may be practiced without these specific details. Thedescriptions and representations herein are the common means used bythose experienced or skilled in the art to most effectively convey thesubstance of their work to others skilled in the art. In otherinstances, well-known methods, procedures, components, and circuitryhave not been described in detail to avoid unnecessarily obscuringaspects of the present invention.

Reference herein to “one embodiment” or “an embodiment” means that aparticular feature, structure, or characteristic described in connectionwith the embodiment can be included in at least one embodiment of theinvention. The appearances of the phrase “in one embodiment” in variousplaces in the specification are not necessarily all referring to thesame embodiment, nor are separate or alternative embodiments mutuallyexclusive of other embodiments. Further, the order of blocks in processflowcharts or diagrams representing one or more embodiments of theinvention do not inherently indicate any particular order nor imply anylimitations in the invention.

Embodiments of the present invention are discussed herein with referenceto FIGS. 1A-9. However, those skilled in the art will readily appreciatethat the detailed description given herein with respect to these figuresis for explanatory purposes as the invention extends beyond theselimited embodiments.

The present invention is directed to methods and systems for numericalsimulation of structural behaviors using ME-FEM. One of the techniquesused in ME-FEM is called generalized meshfree (GMF) approximation orprocedure.

Generalized Meshfree (GMF) Approximation

A first-order GMF convex approximation is introduced to approximate thedisplacement field for ME-FEM. Implementation of meshfree convexapproximation imposes the weak Kronecker-delta property at theboundaries of each finite element and therefore the essential boundarycondition can be easily treated. Fundamental idea of the GMFapproximation is to introduce an enriched basis function in the Shepardfunction and its satisfaction of linear consistency. The selection ofthe enriched basis function determines whether the GMF approximation hasconvexity property. A convex hull conv(Λ) of a node set Λ={x_(i),i=1, .. . n}⊂

² is defined by

$\begin{matrix}{{{conv}(\Lambda)} = \left\{ {\left. {\sum\limits_{i = 1}^{n}\;{\alpha_{i}x_{i}}} \middle| {x_{i} \in \Lambda} \right.,{\alpha_{i} \in {\mathbb{R}}},{\alpha_{i} \geq 0},{{\sum\limits_{i = 1}^{n}\;\alpha_{i}} = 1},{i = 1},2,\ldots} \right\}} & (12)\end{matrix}$The GMF procedure is used for constructing convex approximations of agiven smooth function u, that has the form

${u^{h}(x)} = {\sum\limits_{i = 1}^{n}\;{{\Psi_{i}(x)}u_{i}}}$with the generating function Ψ_(i): conv(Λ)→R satisfying a polynomialreproduction property. An exemplary polynomial reproduction property isexpressed as

${\sum\limits_{i = 1}^{n}\;{{\Psi_{i}(x)}x_{i}}} = {x{\forall{x \in {{{conv}(\Lambda)}.}}}}$First-order GMF approximation in multi-dimension can be written asfollows:

$\begin{matrix}{{\Psi_{i}\left( {x,\lambda_{r}} \right)} = {\frac{\psi_{i}}{\psi} = \frac{{\phi_{a}\left( {x;X_{i}} \right)}{\Gamma_{i}\left( {X_{i},\lambda_{r}} \right)}{B\left( {X_{i},\lambda_{r},\Theta} \right)}}{\sum\limits_{j = 1}^{n}\;{{\phi_{a}\left( {x;X_{j}} \right)}{\Gamma_{j}\left( {X_{j},\lambda_{r}} \right)}{B\left( {X_{j},\lambda_{r},\Theta} \right)}}}}} & (13)\end{matrix}$subject to linear constraints

$\begin{matrix}{{{R_{r}\left( {x,\lambda_{r}} \right)} = {{\sum\limits_{i = 1}^{n}\;{\Psi_{i}X_{i}}} = 0}}{{where}\text{:}}} & (14) \\{\psi_{i} = {{\phi_{a}\left( {x;X_{i}} \right)}{\Gamma_{i}\left( {X_{i},\lambda_{r}} \right)}{B\left( {X_{i},\lambda_{r},\Theta} \right)}}} & (15) \\{\psi = {{\sum\limits_{i = 1}^{n}\;\psi_{i}} = {\sum\limits_{i = 1}^{n}\;{{\phi_{a}\left( {x;X_{i}} \right)}{\Gamma_{i}\left( {X_{i},\lambda_{r}} \right)}{B\left( {X_{i},\lambda_{r},\Theta} \right)}}}}} & (16) \\{X_{i} = {x - x_{i}}} & (17)\end{matrix}$

-   -   φ_(a) (x; X_(i)): weight function of node i with support size        supp(φ_(a)(x; X_(i)))=a_(i),    -   Γ_(i)(X_(i), λ_(r)): the basis function of the GMF        approximation,    -   x: the coordinate of interior point (i.e., a fixed point),    -   x_(i): the coordinate of node i,    -   n: the number of nodes within the support size a(x) at the fixed        point x,    -   λ_(r)(x) (r=1, 2, . . . , m): constraint parameters which have        to be decided,    -   m: the number of constraints (m=1 in one dimension, m=2 in        two-dimension, and m=3 in three-dimension),    -   B(X_(i),λ_(r),Θ): a blending function and Θ is internal        functions or variables.

In the GMF approximation, the property of the partition of unity isautomatically satisfied with the normalization in Eq. (13). Blendingfunction B(X_(i),λ_(r),Θ) is introduced to modify the basis function forthe generation of a particular GMF approximation with special purposes,for example, satisfying the weak Kronecker-delta property for thenon-convex approximation. In constructing the first-order convexapproximation, blending function B(X_(i),λ_(r),Θ) is a unity.Furthermore, the GMF approximation is completed by finding λ to satisfyEq. (14), which is a constraint equation guaranteeing the GMFapproximation to have linear consistency. To determine λ at any fixedpoint x in Eq. (13), a root-finding algorithm is required and used fornon-linear basis functions.

Convexity of the GMF approximation is determined by the selection of apositive basis function in Eq. (13). Table 1 summarizes various basisfunctions typically used for generating either convex or non-convex GMFapproximations. A combination of basis functions or other monotonicallyincreasing or decreasing functions can also be considered as other kindsof GMF basis functions.

The GMF(exp) approximation listed in Table 1 is known as themeshfree-enriched (ME) approximation with Shannon entropy, while theGMF(MLS) approximation is as the Moving Least Squares (MLS) or theReproducing Kernel (RK) approximation. Additionally, the GMF(Renyi)approximation is the ME approximation with Renyi entropy. Moreover, theGMF approximation can be extended to the higher-order, which can be usedfor improving accuracy of meshfree methods. A convex GMF approximationconstructed using an inverse tangent basis function is employed in thedisplacement approximation for the ME-FEM method and is denoted byGMF(atan) in Table 1. A cubic spline window function with a rectangularsupport is chosen as the weight function in Eq. (15). The rectangularsupport is used because it can fit the isoparametric mapping domain of aquadrilateral element.

TABLE 1 Examples of basis functions in the GMF approximations ConvexityBasis function Abbreviation Note Convex e^(x) GMF(exp) ME Approximationapproximation (Shannon entropy) 1 + tanh(x) GMF(tanh) New approximation$1 + {\frac{2}{\pi}{\tan^{- 1}(x)}}$ GMF(atan) New approximationNon-convex approximation$\left( {1 + {\frac{1 - \alpha}{\alpha}x}} \right)^{\frac{1}{\alpha - 1}}\mspace{14mu}\left( {{Renyi}\mspace{14mu}{basis}\mspace{14mu}{function}} \right)$GMF(Renyi) ME Approximation (Renyi entropy) (0.5 < α < 10.5 < α < 1) MLSapproximation (α = 2α = 2) 1 + x GMF(MLS) MLS approximation 1 + x³GMF(x³) New approximation e^(x)(1 + x³) GMF(exp · x³) New approximation

Enriched Finite Element Interpolation by GMF Approximation

Let M_(h)=∪_(e)Q_(e) be the regular division (sometimes referred to astriangulation) of a domain Ω 202 shown in FIG. 2A into convexquadrilateral elements 204 (only few shown). Each element Q_(e) 204contains four corner nodes N_(i), 1≦i≦4 206 a-d and at least onemeshfree enriched (ME) node 208 K_(i), 1≦i≦nb. The quantity nb is thenumber of enriched nodes in the element Q_(e) 204. The element Q_(e) 204can also be transformed into an image of a reference square 214 Q_(e)[−1,1]x[−1,1] through an isoparametric mapping represented in avector form as follows:x=F _(e)(ξ)  (18)where x=[x,y]^(T) and ξ=[ξ,η]^(T).The mapping function F_(e) 220 maps the reference square Q _(e) 214 to ageneral irregular element Q_(e) 204. For illustration and descriptionsimplicity, a five-noded meshfree enriched (ME) finite element 204 withthe fifth node K₁ 208 located at the center of the element is shown inFIG. 2B (i.e., nb is equal to one (1) for this example). Variousprocedures of meshfree enrichment using different number of nodes anddifferent meshfree convex approximation can be facilitated to achievesubstantially similar objective.

Because shape functions of an element are constructed by the meshfreeconvex approximation, the shape functions exhibit the followingconvexity propertiesΨ_(i)(ξ_(j))=δ_(ij)1≦i,j≦4Ψ₅(ξ)<0∀ξε Q _(e)Ψ₅(ξ)=0∀ξε∂ Q _(e)  (19)

Assuming that each node is assigned to a weight function with samerectangular support supp(φ_(a)(ξ;ξ_(i)))=a_(i)=a in the reference squareQ _(e) 214, the tensor product weight function φ_(a)(x;X_(i)) in Eq.(15) is defined as follows:

$\begin{matrix}{{{\phi_{a}\left( {\xi;\xi_{i}} \right)}:\left. {\overset{\_}{Q}}_{e}\rightarrow{\left\lbrack {0,1} \right\rbrack{\phi_{a}\left( {\xi;\xi_{i}} \right)}} \right.} = {{\phi\left( \frac{\xi - \xi_{i}}{a} \right)}{\phi\left( \frac{\eta - \eta_{i}}{a} \right)}}} & (20)\end{matrix}$

Based on the above assumption, there exist certain symmetry propertiesin the shape functions and their derivatives. Using the five-nodedME-FEM shape functions, the isoparametric mapping function F_(e)(ξ) 220in Eq. (18) can be defined as follows:

$\begin{matrix}{\mspace{79mu}{{F_{e}:\left. {\overset{\_}{Q}}_{e}\rightarrow Q_{e} \right.},{x = {{F_{e}(\xi)} = {\left( {F_{e}^{1},F_{e}^{2}} \right) = {{\left( {{\sum\limits_{i = 1}^{5}\;{x_{i}{\Psi_{i}\left( {\xi,\eta} \right)}}},{\sum\limits_{i = 1}^{5}\;{y_{i}{\Psi_{i}\left( {\xi,\eta} \right)}}}} \right)\mspace{14mu}{for}\mspace{14mu}{all}\mspace{14mu}\xi} \in {\overset{\_}{Q}}_{e}}}}}}} & (21)\end{matrix}$

With the above notations, the following approximation space for thedisplacement field is then defined.V ^(h)(Ω)={u ^(h) :u ^(h) ε[H ₀ ¹]² ,u ^(h)|_(Q) _(e) =ū ^(h) ∘F _(e) ⁻¹,ū ^(h) εP ₁( Q _(e))for all Q _(e) εM _(h)}  (22)where P₁( Q _(e))=span {Ψ_(i),i=1 . . . 5} denotes the space contains aset of basis functions in reference square Q _(e) 214.

The fifth or ME node in the ME finite element provides a functionsimilar to the bubble function in the incompatible or enhanced strainelements in the traditional finite element method. However, the basisfunction associated with the fifth node is not orthogonal to the basisfunction associated with the corner node due to its dependence ofmeshfree shape functions. Thus the fifth node cannot be eliminated bythe static condensation procedure in a ME finite element. Unlike theunknown coefficients of the traditional or classical bubble function,the fifth node carries physical displacements. In computing the MEfinite element's body force vectors, the fifth node must be considered.According to Eq. (21), the isoparametric mapping at the fifth node reads

$\begin{matrix}{\left( {x_{5},y_{5}} \right) = {{F_{e}\left( {0,0} \right)} = \left( {{\sum\limits_{i = 1}^{5}\;{x_{i}{\Psi_{i}\left( {0,0} \right)}}},{\sum\limits_{i = 1}^{5}\;{y_{i}{\Psi_{i}\left( {0,0} \right)}}}} \right)}} & (23)\end{matrix}$Consequently, the global co-coordinates of the fifth node is determinedby

$\begin{matrix}{\left( {x_{5},y_{5}} \right) = {\left( {\frac{\sum\limits_{i = 1}^{4}\;{{\Psi_{i}\left( {0,0} \right)}x_{i}}}{\left( {1 - {\Psi_{5}\left( {0,0} \right)}} \right)},\frac{\sum\limits_{i = 1}^{4}\;{{\Psi_{i}\left( {0,0} \right)}y_{i}}}{\left( {1 - {\Psi_{5}\left( {0,0} \right)}} \right)}} \right) = {\left( {\frac{\sum\limits_{i = 1}^{4}\;{{\Psi_{i}\left( {0,0} \right)}x_{i}}}{\sum\limits_{i = 1}^{4}\;{\Psi_{i}\left( {0,0} \right)}},\frac{\sum\limits_{i = 1}^{4}\;{{\Psi_{i}\left( {0,0} \right)}y_{i}}}{\sum\limits_{i = 1}^{4}\;{\Psi_{i}\left( {0,0} \right)}}} \right) = \left( {\frac{\sum\limits_{i = 1}^{4}\; x_{i}}{4},\frac{\sum\limits_{i = 1}^{4}\; y_{i}}{4}} \right)}}} & (24)\end{matrix}$which is the centroid of the quadrilateral (i.e., the 5-node ME finiteelement 204).

For a 4-noded meshfree-enriched finite element (triangular element 224having three corner nodes shown as solid dots and one ME node shown as acircle), an exemplary isoparametric mapping scheme is shown in FIG. 2C.It can be seen circular supports 230 are used for the triangular element224 instead of the rectangular support for the quadrilateral element204.

Further, the present invention can apply to a three-dimensional element(i.e., tetrahedral element 244) shown in FIG. 2D, in which an exemplaryisoparametric mapping scheme is illustrated.

FIG. 3A shows an exemplary shape function of one of the corner nodes 206a-d, while FIG. 3B shows an exemplary shape function of the meshfreeenriched (ME) node 208 (i.e., middle node or fifth node). These shapefunctions are results of GMF(atan) approximation in a five-noded ME-FEMelement 204 which shows the Kronecker-delta property at the boundary.The shape functions of the ME-FEM element reduce to the standardbilinear finite element shape functions along the element edge thatexhibit the convexity approximation property. It is noted that thederivatives of the shape functions (listed in Eq, (79) below) in thefive-noded ME-FEM element 204 do not degenerate to the finite elementderivatives along the element edge due to the non-zero value ofderivatives of the ME or fifth node 208.

The Jacobian matrix J_(ξ) is defined to be the gradient of the mappingfunction F_(e) by

$\begin{matrix}{{J_{ij} = \frac{\partial{F_{e}^{i}(\xi)}}{\partial\xi_{j\;}}},i,{j = 1},2} & (25)\end{matrix}$or in matrix formJ _(ξ) =G _(F) G _(F) ^(T)  (26)where matrices G_(F) and G_(x) are defined as follows:

$\begin{matrix}{G_{F} = {{\nabla_{({\xi,\eta})}\left\lbrack {\Psi_{i}\left( {\xi,\eta} \right)} \right\rbrack_{e}} = \begin{bmatrix}\frac{\partial\Psi_{1}}{\partial\xi} & \ldots & \frac{\partial\Psi_{5}}{\partial\xi} \\\frac{\partial\Psi_{1}}{\partial\eta} & \ldots & \frac{\partial\Psi_{5}}{\partial\eta}\end{bmatrix}}} & (27) \\{G_{x} = \begin{bmatrix}x_{1} & \ldots & x_{5} \\y_{1} & \ldots & y_{5}\end{bmatrix}} & (28) \\{\nabla_{({\xi,\eta})}{= \begin{bmatrix}\frac{\partial}{\partial\xi} \\\frac{\partial}{\partial\eta}\end{bmatrix}}} & (29)\end{matrix}$

Analytical proof of invertibility of the mapping function F_(e) (ornon-singularity of the Jacobian) is not an easy task due to the rationalnature of the meshfree shape functions and should be addressed by theinverse function theorem. In engineering applications, it is more commonto examine the invertibility at the Gauss integration points used forquadrature of the ME element 204.

The isoparametric mapping in the five-noded ME-FEM quadrilateral element204 preserves the global element area A_(e) numerically, where

$A_{e} = {{\int_{\Omega_{e}}{\mathbb{d}\Omega}} = {{\int_{- 1}^{1}{\int_{- 1}^{1}{{\det(J)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}} = {\overset{4}{\sum\limits_{J = 1}}{{\det\left( J_{\xi\;} \right)}_{J}.}}}}$In other words, the sum of det(J_(ξ)) using a 2×2 Gaussian integrationscheme preserves the area of global element Q_(e) 204 despite the ME-FEMelement shape functions and their derivatives are rational functions andthe nodal support sizes are adjustable.

Proof

From Eq. (26), the area of global element Q_(e) is obtained using theconventional 2×2 Gaussian integration given by

$\begin{matrix}\begin{matrix}{A_{e} = {\sum\limits_{J = 1}^{4}{\det\left( J_{\xi} \right)}_{J}}} \\{= {\sum\limits_{J = 1}^{4}{\det\left( {{G_{F}\left( \xi_{gJ} \right)}G_{x}^{T}} \right)}}} \\{= {\sum\limits_{J = 1}^{4}\begin{pmatrix}{{\sum\limits_{I = 1}^{5}{\frac{\partial{\Psi_{I}\left( \xi_{gJ} \right)}}{\partial\xi}{x_{I} \cdot {\sum\limits_{I = 1}^{5}{\frac{\partial{\Psi_{I}\left( \xi_{gJ} \right)}}{\partial\eta}y_{I}}}}}} -} \\{\sum\limits_{I = 1}^{5}{\frac{\partial{\Psi_{I}\left( \xi_{gJ} \right)}}{\partial\xi}{y_{I} \cdot {\sum\limits_{I = 1}^{5}{\frac{\partial{\Psi_{I}\left( \xi_{gJ} \right)}}{\partial\eta}x_{I}}}}}}\end{pmatrix}}}\end{matrix} & (30)\end{matrix}$The first-order consistency condition gives

$\begin{matrix}{{{{G_{F}\left( \xi_{gi} \right)}G_{\xi}^{T}} = {{I\mspace{14mu}{for}\mspace{14mu} i} = 1}},2,3,4} & (31) \\{G_{\xi} = \begin{bmatrix}{- 1} & 1 & {- 1} & {- 1} & 0 \\{- 1} & {- 1} & 1 & 1 & 0\end{bmatrix}} & (32)\end{matrix}$Using the symmetry properties in Eq. (79) together with the first-orderconsistency conditions leads to the following explicit expression ofA_(e)

$\begin{matrix}{A_{e} = \frac{\begin{matrix}{{\left( {{- x_{1}} + x_{2} + x_{3} - x_{4}} \right) \cdot \left( {{- y_{1}} - y_{2} + y_{3} + y_{4}} \right)} -} \\{\left( {{- y_{1}} + y_{2} + y_{3} - y_{4}} \right) \cdot \left( {{- x_{1}} - x_{2} + x_{3} + x_{4}} \right)}\end{matrix}}{2}} & (33)\end{matrix}$A_(e) is the exact area of global element Q_(e) 204 and is independentof co-coordinates or location of the ME or fifth node 208 and nodalsupport size a. Furthermore, the ME-FEM formulation preserves the linearexactness in the Galerkin approximation of Dirichlet boundary valueproblem. i.e.∫_(Ω) B _(I) dΩ=0 or ∫_(Ω) b _(I) dΩ=0 ∀ interior node IεΩ∪Γ  (34)which is known as integration constraints (IC) in meshfree Galerkinmethod. In Eq. (34), B_(I) is the gradient matrix which is given by

$\begin{matrix}{{{B_{I}(\xi)} = \begin{bmatrix}\frac{\partial{\Psi_{I}(\xi)}}{\partial x} & 0 \\0 & \frac{\partial{\Psi_{I}(\xi)}}{\partial y} \\\frac{\partial{\Psi_{I}(\xi)}}{\partial y} & \frac{\partial{\Psi_{I}(\xi)}}{\partial x}\end{bmatrix}}{and}{{b_{I}(\xi)} = \begin{bmatrix}\frac{\partial{\Psi_{I}(\xi)}}{\partial x} \\\frac{\partial{\Psi_{I}(\xi)}}{\partial y}\end{bmatrix}}} & (35)\end{matrix}$is a vector of the shape function derivatives of node I.

Proof

First consider the case that interior node is the center node K₁ 208 ofthe global element Q_(e) 204. In each element,

$\begin{matrix}{{(36){\begin{matrix}{\;{{\int_{Q_{e}}{b_{5}{\mathbb{d}\Omega}}} = {\int_{- 1}^{1}{\int_{- 1}^{1}{b_{5}{\det\left( J_{\xi} \right)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}}} \\{= {\int_{- 1}^{1}{\int_{- 1}^{1}{{J_{\xi}^{- 1} \cdot {\overset{\_}{b}}_{5}}{\det\left( J_{\xi} \right)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}} \\{= {\int_{- 1}^{1}{\int_{- 1}^{1}{{J_{\xi}^{*} \cdot {\overset{\_}{b}}_{5}}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}}\end{matrix}\mspace{110mu} = {\begin{bmatrix}{\sum\limits_{j = 1}^{4}\begin{pmatrix}{\;{{\left( {\sum\limits_{i = 1}^{5}\;{y_{i}\;\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\eta}}} \right)\frac{\partial{\Psi_{5}\left( \xi_{gj} \right)}}{\partial\xi}} -}} \\{\left( {\sum\limits_{i = 1}^{5}\;{y_{i}\frac{\partial{\Psi_{i}\left( \xi_{gi} \right)}}{\partial\xi}}} \right)\frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\eta}}\end{pmatrix}} \\{\sum\limits_{j = 1}^{4}\begin{pmatrix}{\;{{\left( {\sum\limits_{i = 1}^{5}\;{x_{i}\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\eta}}} \right)\frac{\partial{\Psi_{5}\left( \xi_{gj} \right)}}{\partial\xi}} -}} \\{\left( {\sum\limits_{i = 1}^{5}\;{x_{i}\frac{\partial{\Psi_{i}\left( \xi_{gi} \right)}}{\partial\xi}}} \right)\frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\eta}}\end{pmatrix}}\end{bmatrix}\mspace{110mu} = {\begin{bmatrix}{\sum\limits_{i = 1}^{5}\left( {y_{i}{\sum\limits_{j = 1}^{4}\;\left( \;\begin{matrix}{{\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\eta}\frac{\partial{\Psi_{5}\left( \xi_{gj} \right)}}{\partial\xi}} -} \\{\frac{\partial{\Psi_{i}\left( \xi_{gi} \right)}}{\partial\xi}\frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\eta}}\end{matrix} \right)}}\; \right)} \\{\;{\sum\limits_{i = 1}^{5}\left( {x_{i}{\sum\limits_{j = 1}^{4}\;\begin{pmatrix}{{\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\eta}\frac{\partial{\Psi_{5}\left( \xi_{gj} \right)}}{\partial\xi}} -} \\{\frac{\partial{\Psi_{i}\left( \xi_{gi} \right)}}{\partial\xi}\frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\eta}}\end{pmatrix}}} \right)}}\end{bmatrix}\mspace{110mu} = \left\lbrack \begin{matrix}{\sum\limits_{i = 1}^{5}\;{y_{i}g_{i5}}} \\{\sum\limits_{i = 1}^{5}\;{x_{i}g_{i5}}}\end{matrix} \right\rbrack}}}\mspace{11mu}{where}\mspace{709mu}(37){\overset{\_}{b_{I}} = {{\nabla_{({\xi,\eta})}\Psi_{I}} = \begin{bmatrix}\frac{\partial{\Psi_{I}\left( {\xi,\eta} \right)}}{\partial\xi} \\\frac{\partial{\Psi_{I}\left( {\xi,\eta} \right)}}{\partial\eta}\end{bmatrix}}}\mspace{706mu}(38){{J_{\xi}^{*}(\xi)} = {\det\left( J_{\xi} \right)}}\;\mspace{706mu}(39){g_{iI} = {\sum\limits_{j = 1}^{4}\;\begin{pmatrix}{{\frac{\partial{\Psi_{i}\left( \xi_{gJ} \right)}}{\partial\eta}\frac{\partial{\Psi_{I}\left( \xi_{gj} \right)}}{\partial\xi}} -} \\{\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\xi}\frac{\partial{\Psi_{I}\left( \xi_{gj} \right)}}{\partial\eta}}\end{pmatrix}}}}} & \;\end{matrix}$Using the symmetry properties in Eq. (79), a calculation leads to∫_(Q) _(e) b ₅ dΩ=0  (40)Let's next consider an interior node I=N₃ ^(e1) connected by severalneighboring elements where the superscript “e1” denotes the elementnumber and subscript “3” is the node number of that particular elementshown in FIG. 4. Similar to the above derivation, we have

$\begin{matrix}\begin{matrix}{{\int_{Q_{e}}{b_{3\;}^{e\; 1}{\mathbb{d}\Omega}}} = {\int_{- 1}^{1}{\int_{- 1}^{1}{b_{3}^{e\; 1}{\det\left( J_{\xi}^{e\; 1} \right)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}} \\{= \begin{bmatrix}{\sum\limits_{j = 1}^{4}\begin{pmatrix}{{\left( {\sum\limits_{i = 1}^{5}{y_{i}^{e\; 1}\frac{\partial{\Psi_{i}\left( \xi_{{gj}\;} \right)}}{\partial\eta}}} \right)\frac{\partial{\Psi_{3}\left( \xi_{gj} \right)}}{\partial\xi}} -} \\{\left( {\sum\limits_{i = 1}^{5}{y_{i}^{e\; 1}\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\xi}}} \right)\frac{\partial{\Psi_{3}\left( \xi_{gj} \right)}}{\partial\eta}}\end{pmatrix}} \\{\sum\limits_{j = 1}^{4}\begin{pmatrix}{{\left( {\sum\limits_{i = 1}^{5}{x_{n}^{e\; 1}\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\eta}}} \right)\frac{\partial{\Psi_{3}\left( \xi_{gj} \right)}}{\partial\xi}} -} \\{\left( {\sum\limits_{i = 1}^{5}{x_{i}^{e\; 1}\frac{\partial{\Psi_{i}\left( \xi_{gj} \right)}}{\partial\xi}}} \right)\frac{\partial{\Psi_{3}\left( \xi_{gj} \right)}}{\partial\eta}}\end{pmatrix}}\end{bmatrix}} \\{= \begin{bmatrix}{{y_{2}^{e\; 1}g_{23}^{e\; 1}} + {y_{4}^{e\; 1}g_{43}^{e\; 1}} + {y_{5}^{e\; 1}g_{53}^{e\; 1}}} \\{{x_{2}^{e\; 1}g_{23}^{e\; 1}} + {x_{4}^{e\; 1}g_{43}^{e\; 1}} + {x_{5}^{e\; 1}g_{53}^{e\; 1}}}\end{bmatrix}}\end{matrix} & (41)\end{matrix}$The result in Eq. (41) is independent of the co-coordinates of node N₃^(e1) and its diagonal node N₁ ^(e1) in element e1. For adjacentelements e1 and e2,

$\begin{matrix}{{{\int_{Q_{e\; 1}}{b_{3}^{e\; 1}{\mathbb{d}\Omega}}} + {\int_{Q_{e\; 2}}{b_{4}^{e\; 2}{\mathbb{d}\Omega}}}} = {\quad{\begin{bmatrix}\begin{matrix}{{y_{2}^{e\; 1}g_{23}^{e\; 1}} + {y_{4}^{e\; 1}g_{43}^{e\; 1}} +} \\{y_{5}^{e\; 1}g_{53}^{e\; 1}}\end{matrix} \\\begin{matrix}{{x_{2}^{e\; 1}g_{23}^{e\; 1}} + {x_{4}^{e\; 1}g_{43}^{e\; 1}} +} \\{x_{5}^{e\; 1}g_{53}^{e\; 1}}\end{matrix}\end{bmatrix} + \begin{bmatrix}\begin{matrix}{{y_{1}^{e\; 2}g_{14}^{e\; 2}} + {y_{3}^{e\; 2}g_{34}^{e\; 2}} +} \\{y_{5}^{e\; 2}g_{54}^{e\; 2}}\end{matrix} \\\begin{matrix}{{x_{1}^{e\; 2}g_{14}^{e\; 2}} + {x_{3}^{e\; 2}g_{34}^{e\; 2}} +} \\{x_{5}^{e\; 2}g_{54}^{e\; 2}}\end{matrix}\end{bmatrix}}}} & (42)\end{matrix}$Using the relationships x₂ ^(e1)=x₁ ^(e2), y₂ ^(e1)=y₁ ^(e2) and g₂₃^(e1)=g₁₄ ^(e2), Eq. (42) is rewritten as

$\begin{matrix}{{{\int_{Q_{e\; 1} + Q_{e\; 2}}b_{I}^{e\; 1}} + {b_{I}^{e\; 2}{\mathbb{d}\Omega}}} = {\begin{bmatrix}{{y_{4}^{e\; 1}g_{43}^{e\; 1}} + {y_{5}^{e\; 1}g_{53}^{e\; 1}}} \\{{x_{4}^{e\; 1}g_{43}^{e\; 1}} + {x_{5}^{e\; 1}g_{53}^{e\; 1}}}\end{bmatrix} + \begin{bmatrix}{{y_{3}^{e\; 2}g_{34}^{e\; 2}} + {y_{5}^{e\; 2}g_{54}^{e\; 2}}} \\{{x_{3}^{e\; 2}g_{34}^{e\; 2}} + {x_{5}^{\;{e\; 2}}g_{54}^{e\; 2}}}\end{bmatrix}}} & (43)\end{matrix}$It is noted that I=N₃ ^(e1)=N₄ ^(e2) from FIG. 4. Subsequently, formultiple elements n1 one can easily verify that

$\begin{matrix}{{\int_{\Omega}{b_{I}{\mathbb{d}\Omega}}} = {{\sum\limits_{n = 1}^{nl}{\int_{- 1}^{1}{\int_{- 1}^{1}{b_{I}^{n}{\det\left( J_{\xi}^{n} \right)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}} = 0}} & (44)\end{matrix}$

Locking-Free Meshfree-Enriched Finite Element FormulationDivergence-Free Meshfree-Enriched Finite Element Interpolation

From the nearly incompressible problem in Eq. (9), the error between theexact solution u and the unique solution u^(h) is obtained from astandard displacement-based Galerkin method is estimated and expressedas the energy error norm by

$\begin{matrix}{{{{u - u^{h}}}_{E}^{2} \leq {C_{0}\inf\limits_{v^{h} \in V^{h}}{{u - v^{h}}}_{E}^{2}} \leq {C_{0}{{u - {\prod_{h}u}}}_{E}^{2}}} = {{C_{0}{A\left( {{u - {\prod_{h}u}},{u - {\prod_{h}u}}} \right)}} = {C_{0}\left( {{2\mu{{u - {\prod_{h}u}}}^{2}} + {\lambda{{\nabla{\cdot \left( {u - {\prod_{h}u}} \right)}}}_{0,\Omega}^{2}}} \right)}}} & (45)\end{matrix}$The energy norm in Eq. (45) is defined by

$\begin{matrix}{{u}_{E} = \left( {A\left( {u,u} \right)} \right)^{\frac{1}{2}}} & (46)\end{matrix}$Π_(h):H¹(Ω)→V^(h)(Ω) is an interpolation operator associated with theisoparametric mapping F_(e) 220 in Eq. (21) and defined byΠ_(h) u= Π _(h)(u∘F _(e))∘F _(e) ⁻¹  (47)where Π _(h):H¹( Q _(e))→P₁( Q _(e)) is a projection operator thatpreserves polynomials of degree≦1; i.e., Π _(h)p=p, ∀pεP₁( Q _(e)) andP₁( Q _(e)) is the space of polynomials of degree≦1 on Q _(e). ∥•∥_(m)and |•|_(m) are Sobolev norm and semi-norm of order m, respectively.Constant C₀ in Eq. (45) does not dependent on the element size h and thefunctions considered. It is noted that the well-known Céa's lemma forthe first inequality in Eq. (45) is used to estimate the error. Let'sfurther assume the regularity on the solution uεH²(Ω) and use theapproximation property for the first-order approximation to yield theglobal approximation error

$\begin{matrix}{{{u - u^{h}}}_{E}^{2} \leq {C_{0}\left( {{2\mu{{u - {\prod_{h}u}}}_{1,\Omega}^{2}} + {\lambda{{\nabla{\cdot \left( {u - {\prod_{h}u}} \right)}}}_{0,\Omega}^{2}}} \right)} \leq {{C_{01}\mu\; h^{2}{u}_{2,\Omega}^{2}} + {C_{10}\lambda\; h^{2}{{\nabla{\cdot u}}}_{1,\Omega}^{2}}} \leq {C_{11}{h^{2}\left( {{u}_{2,\Omega}^{2} + {\lambda{{\nabla{\cdot u}}}_{1,\Omega}^{2}}} \right)}}} & (48)\end{matrix}$

From Eq. (48) it is clear that solution u^(h) will converge to the exactsolution u as mesh is refined in the compressible case. However theconvergence of u^(h) may not be achieved by a low-order approximation asλ→∞. For the energy error norm to be bounded in the incompressiblelimit, the divergence-free condition has to be enforced strongly using ahigh-order approximation or enforced weakly by a mixed/reducedintegration/modified variational form strategy and guarded by theinf-sup condition.

Instead of enforcing the divergence-free condition weakly, thedivergence-free condition is enforced strongly but point-wise at theGaussian quadrature points 410 (shown as triangles in FIG. 4). That is,∇·u ^(h)(ξ_(gi))|_(Q) _(e) →0 as λ→∞∀Q _(e) εM _(h) ,i=1,2,3,4 are Gausspoints  (49)

The approximation of the divergence ∇·u^(h) is evaluated at the Gausspoints 410, which can be expressed by∇·u ^(h)=(ξ_(gi))=tr(∇u ^(h)(ξ_(gi)))=tr(J _(ξ) ⁻¹·∇_((ξ,η)) u^(h)(ξ_(gi)))=tr((G _(F) G _(x) ^(T))⁻¹·∇_((ξ,η)) u ^(h)(ξ_(i)))  (50)where ξ_(gi)=(ξ_(gi),η_(gi)), i=1, 2, 3, 4 are Gauss points 410. Theterm ∇_((ξ,η))u^(h)(ξ_(gi)) in Eq. (50) is rewritten as

$\begin{matrix}{{{\nabla_{({\xi,\eta})}{u^{h}\left( \xi_{gi} \right)}} = {{\begin{bmatrix}\frac{\partial{\Psi_{1}\left( \xi_{gi} \right)}}{\partial\xi} & \ldots & \frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\xi} \\\frac{\partial{\Psi_{1}\left( \xi_{gi} \right)}}{\partial\eta} & \ldots & \frac{\partial{\Psi_{5}\left( \xi_{gi} \right)}}{\partial\eta}\end{bmatrix}\begin{bmatrix}u_{1} & v_{1} \\\vdots & \vdots \\u_{5} & v_{5}\end{bmatrix}} = {{G_{F}\left( \xi_{gi} \right)}u_{e}}}},\mspace{20mu}{i = 1},2,3,4} & (51)\end{matrix}$where u_(e) is a collection of nodal displacement for the element Q_(e)204 given by

$\begin{matrix}{u_{e}^{T} = \begin{bmatrix}u_{1} & \ldots & u_{5} \\v_{1} & \ldots & v_{5}\end{bmatrix}} & (52)\end{matrix}$

To show that the divergence-free condition can be achieved point-wiseusing the ME-FEM approximation, the eigenvalue analysis is used. Using5-noded element 204 Q_(e)=[−1,1]×[−1,1], the global co-coordinatescoincide with the reference co-coordinates. In other words, the Jacobianmatrix J_(ξ) in Eq. (26) becomes J_(ξ)=G_(F)G_(x) ^(r)=I and theexpression of divergence ∇·u^(h) can be degenerated to a simple form.FIG. 5A illustrates one of the volumetric locking modes 510 in arectangular bilinear four node finite element where the displacementfields are given in the global co-coordinates by[u,v]=[c ₁ xy,0]  (53)

In Eq. (53), c₁ is a coefficient describing the deformation of thelocking mode which is a function of nodal displacements. The divergence∇·u^(h) is not null evaluating at all four Gauss points 410 andinterpolation does not meet the divergence-free condition of Eq. (49) inthe incompressible limit. On the other hand, the divergence ∇·u^(h) insame deformation mode using the ME-FEM approximation evaluating at Gausspoints 410 can be expressed using Eq. (50) as

$\begin{matrix}\begin{matrix}{{\nabla{\cdot {u^{h}\left( \xi_{gi} \right)}}} = {{tr}\left( {J_{\xi}^{- 1} \cdot {\nabla_{({\xi,\eta})}{u^{h}\left( \xi_{gi} \right)}}} \right)}} \\{= {{tr}\left( {I \cdot {\nabla_{({\xi,\eta})}{u^{h}\left( \xi_{gi} \right)}}} \right)}} \\{= {{\sum\limits_{k = 1}^{5}{{\Psi_{{k,\xi}\;}\left( {\xi_{gi},\eta_{gi}} \right)}u_{k}}} +}} \\{\sum\limits_{k = 1}^{5}{{\Psi_{k,\eta}\left( {\xi_{gi},\eta_{gi}} \right)}v_{k}}}\end{matrix} & (54) \\{{{{for}\mspace{14mu} i} = 1},2,3,4} & \;\end{matrix}$

According to Eq. (49), the divergence-free interpolation in ME-FEMelement Q_(e) 204 is∇·u ^(h)(ξ_(gi))=0 for i=1,2,3,4  (55)or in matrix form

$\begin{matrix}{{\begin{bmatrix}{\Psi_{5,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)} & {\Psi_{5,\eta}\;\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)} \\{\Psi_{5,\xi}\left( {\xi_{g\; 2},\eta_{g\; 2}} \right)} & {\Psi_{5,\eta}\left( {\xi_{g\; 2},\eta_{g\; 2}} \right)} \\{\Psi_{5,\xi}\left( {\xi_{g\; 3},\eta_{g\; 3}} \right)} & {\Psi_{5,\eta}\left( {\xi_{g\; 3},\eta_{g\; 3}} \right)} \\{\Psi_{5,\xi}\left( {\xi_{g\; 4},\eta_{g\; 4}} \right)} & {\Psi_{5,\eta}\left( {\xi_{g\; 4},\eta_{g\; 4}}\; \right)}\end{bmatrix}\begin{bmatrix}u_{5} \\v_{5}\end{bmatrix}} = {- \begin{bmatrix}{\sum\limits_{k = 1}^{4}\left( {{{\Psi_{k,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)}u_{k}} + {{\Psi_{k,\eta}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)}v_{k}}} \right)} \\{\sum\limits_{k = 1}^{4}\left( {{{\Psi_{k,\xi}\left( {\xi_{g\; 2},\eta_{g\; 2}} \right)}u_{k}} + {{\Psi_{k,\eta}\left( {\xi_{g\; 2},\eta_{g\; 2}} \right)}v_{k}}} \right)} \\{\sum\limits_{k = 1}^{4}\left( {{{\Psi_{k,\xi}\left( {\xi_{g\; 3},\eta_{g\; 3}} \right)}u_{k}} + {{\Psi_{k,\eta}\left( {\xi_{g\; 3},\eta_{g\; 3}} \right)}v_{k}}} \right)} \\{\sum\limits_{k = 1}^{4}\left( {{{\Psi_{k,\xi}\left( {\xi_{g\; 4},\eta_{g\; 4}} \right)}u_{k}} + {{\Psi_{k,\eta}\left( {\xi_{g\; 4},\eta_{g\; 4}} \right)}v_{k}}} \right)}\end{bmatrix}}} & (56)\end{matrix}$

The above linear system consists of two unknowns (u₅ and v₅) and fourconstraint equations along with the known shape function derivativesevaluated at four Gauss points 410 and eight prescribed nodaldisplacements. However, the above four constraint equations are nottotally independent. Using Eq. (79), we can reduce four constraintequations into two equations and obtain the two unknowns by

$\begin{matrix}{\begin{bmatrix}u_{5} \\v_{5}\end{bmatrix} = \begin{bmatrix}\frac{c_{1}\begin{pmatrix}{{- {\Psi_{1,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)}} - {\Psi_{3,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)} +} \\{{\Psi_{2,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)} + {\Psi_{4,\xi}\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)}}\end{pmatrix}}{\Psi_{5,\eta}\left( {\xi_{g\; 1},\eta_{g\; 1}}\; \right)} \\0\end{bmatrix}} & (57)\end{matrix}$

The above result demonstrates that the enrichment of the ME node inME-FEM element can exactly produce two non-locking modes, which aremissing in the standard bilinear quadrilateral element. As a result, theoverall ME-FEM element behavior is free of volumetric locking.Furthermore, the displacement field of the five-noded ME-FEM element isdivergence-free point-wise at the Gauss points 410. Since the shapefunction derivative Ψ_(5,η)(ξ_(g1),η_(g1)) in Eq. (57) is non-zero(i.e., non-zero everywhere in the element expect at the ME node itself)using the convex approximation, the displacement of the ME node can beuniquely determined. FIG. 5B presents an eigenvalue analysis resultusing the ME-FEM method. The mode shape 520 is free of volumetriclocking A bounded engienvalue near the incompressible limit indicatesthe ME-FEM element is volumetric locking-free. Results from aneigenanalysis also assure that the ME-FEM element contains no spuriouszero-energy mode. Since the volumetric energy term in Eq. (45) isbounded in the incompressible limit, the energy norm error of ME-FEMsolution is expected to be close to the theoretical asymptotic rate ofO(h). Accordingly, the pressure is calculated using the followingformula:p= ^(h) =λ∇·u ^(h)  (58)which remains finite in the incompressible limit, and the L₂—norm errorestimate holds∥∇·(u−u ^(h))∥_(0,Ω) ≦c ₂ h|∇·u| _(1,Ω)  (59)

Although the ME-FEM element does not contain any spurious zero-energymode, the ME-FEM formulation using conventional Gaussian quadrature ruleexperiences pressure oscillation in the near-incompressible analysis forcertain cases. The cause of pressure oscillation in ME-FEM formulationusing Gaussian quadrature rule is different from the one in thediscontinuous-pressure/displacement mixed finite element method whichexhibits a rank-deficiency in the assembled pressure equations.

The pressure oscillation in ME-FEM method is due to the nature of theweighted residual-based Galerkin method using a low-order approximation.Further, the displacement derivatives in ME-FEM method are discontinuousacross element boundaries. A strain smoothing procedure is employed toremove such oscillation. Several post-processing pressure smoothingtechniques can be applied to produce a smooth pressure distribution.However additional computational effort is required and robustperformance in the non-linear analysis may not be guaranteed.

To suppress the possible pressure oscillation, an area-weightedintegration incorporated with the divergence-free ME-FEM interpolationis developed to provide the smoothing effects for strains and pressureas well as divergence-free property.

Area-Weighted Integration

An area-weighted integration based on the strain smoothing concept isintroduced to the ME-FEM element 204 to provide the smoothing effect onstrains and to suppress the pressure oscillation in thenear-incompressible analysis. Additionally, the area-weightedintegration provides divergence-free ME-FEM interpolation in triangularand tetrahedral elements. The strain smoothing is defined as follows:

$\begin{matrix}{{\overset{\sim}{ɛ}}^{h} = {{\overset{\sim}{\nabla}u^{h}} = {\frac{1}{A_{m}}{\int_{\Omega_{m}}{{\nabla u^{h}}{\mathbb{d}\Omega}}}}}} & (60)\end{matrix}$where {tilde over (∇)} denotes the smoothed gradient operator. A_(m) 601is the area of the smoothed domain Ω_(m) 602. Each element 604 containsfour smoothed domains 602 and each smoothed domain 602 has straincontribution from two Gauss points 610 (shown as triangles) shown inFIG. 6A. Accordingly the smoothed strain in the first smoothed domain602 is expressed by

$\begin{matrix}{{\overset{\sim}{ɛ}}^{h} = {\frac{1}{2A_{m}}\left( {{{ɛ\left( {\xi_{g\; 1},\eta_{g\; 1}} \right)}{\det\left( J_{1} \right)}} + {{ɛ\left( {\xi_{g\; 2},\eta_{g\; 2}} \right)}{\det\left( J_{2} \right)}}} \right)}} & (61) \\{A_{m} = {\frac{{\det\left( J_{1} \right)} + {\det\left( J_{2} \right)}}{2} = \frac{{\det\left( {G_{F}G_{x}^{T}} \right)}_{1} + {\det\left( {G_{F}G_{x}^{T}} \right)}_{2}}{2}}} & (62)\end{matrix}$This leads to the smoothed divergence {tilde over (∇)}·u^(h) inone-element model given by

$\begin{matrix}\begin{matrix}{{\overset{\sim}{\nabla}{\cdot u^{h}}} = {{tr}\left( {\overset{\sim}{ɛ}}^{h} \right)}} \\{= {{tr}\left( {\frac{1}{A_{m}}{\int_{\Omega_{m}}{{\nabla u^{h}}{\mathbb{d}\Omega}}}} \right)}} \\{= {\frac{1}{A_{m}}{\int_{\Omega_{m}}{{{tr}\left( {\nabla u^{h}} \right)}{\mathbb{d}\Omega}}}}} \\{= {\frac{1}{A_{m}}{\int_{\Omega_{m}}{{\nabla{\cdot u^{h}}}{\mathbb{d}\Omega}}}}} \\{= \frac{\begin{matrix}{{{\nabla{\cdot {u^{h}\left( \xi_{g\; 1} \right)}}}{\det\left( J_{1} \right)}} +} \\{{\nabla{\cdot {u^{h}\left( \xi_{g\; 2} \right)}}}{\det\left( J_{2} \right)}}\end{matrix}}{{\det\left( J_{1} \right)} + {\det\left( J_{2} \right)}}} \\{= 0}\end{matrix} & (63)\end{matrix}$which verifies the divergence-free condition of Eq. (49) in theincompressible limit.

For adjacent elements 614 a-b shown in FIG. 6B, the smoothing area A_(m)611 is calculated from the sum of four Jacobian determinants in twoelements 614 a-b given by

$\begin{matrix}\begin{matrix}{A_{m} = \frac{{\det\left( J_{2}^{e\; 1} \right)} + {\det\left( J_{3}^{e\; 1} \right)} + {\det\left( J_{1}^{e\; 2} \right)} + {\det\left( J_{4}^{e\; 2} \right)}}{2}} \\{= \frac{\begin{matrix}{{\det\left( {G_{F}G_{x}^{T}} \right)}_{2}^{e\; 1} + {\det\left( {G_{F}G_{x}^{T}} \right)}_{3}^{e\; 1} +} \\{{\det\left( {G_{F}G_{x}^{T}} \right)}_{1}^{e\; 2} + {\det\left( {G_{F}G_{x}^{T}} \right)}_{4}^{e\; 2}}\end{matrix}}{2}}\end{matrix} & (64)\end{matrix}$where superscripts e_(e1) and e₂ are adjacent elements 614 a-b thatshare the same element edge m and two connecting nodes I_(s1) and I_(s2)616 a-b. Similar to the derivation in Eq. (63), it is proved that thesmoothed divergence {tilde over (∇)}·u^(h) in the adjacent elements 614a-b satisfies the divergence-free condition in the incompressible limit.With Eq. (64), the weighted average of the strain-displacement matrix inthe smoothed domain Ω_(m) is given by

$\begin{matrix}{{\overset{\sim}{B}}_{I} = \begin{bmatrix}{\overset{\sim}{b}}_{Ix} & 0 \\0 & {\overset{\sim}{b}}_{Iy} \\{\overset{\sim}{b}}_{Iy} & {\overset{\sim}{b}}_{Ix}\end{bmatrix}} & (65) \\{{\overset{\sim}{b}}_{Ix} = \frac{\begin{matrix}{{\left( b_{Ix} \right)_{2}^{e\; 1}{\det\left( J_{2}^{e\; 1} \right)}} + {\left( b_{Ix} \right)_{3}^{e\; 1}{\det\left( J_{3}^{e\; 1} \right)}} +} \\{{\left( b_{Ix} \right)_{1}^{e\; 2}{\det\left( J_{1}^{e\; 2} \right)}} + {\left( b_{Ix} \right)_{4}^{e\; 2}{\det\left( J_{4}^{e\; 2} \right)}}}\end{matrix}}{2A_{m}}} & (66) \\{{\overset{\sim}{b}}_{Iy} = \frac{\begin{matrix}{{\left( b_{Iy} \right)_{2}^{e\; 1}{\det\left( J_{2}^{e\; 1} \right)}} + {\left( b_{Iy} \right)_{3}^{e\; 1}{\det\left( J_{3}^{e\; 1} \right)}} +} \\{{\left( b_{I\; y} \right)_{1}^{e\; 2}{\det\left( J_{1}^{e\; 2} \right)}} + {\left( b_{y\; x} \right)_{4\;}^{e\; 2}{\det\left( J_{4}^{e\; 2} \right)}}}\end{matrix}}{2A_{m}}} & (67)\end{matrix}$where the components in {tilde over (b)}_(Ix) and {tilde over (b)}_(Iy)are expressed by

$\begin{matrix}{{\begin{Bmatrix}\left( b_{Ix} \right)_{i}^{e\; 1} \\\left( b_{Ix} \right)_{i}^{e\; 1}\end{Bmatrix} = {J_{i}^{e\; 1^{- 1}}\begin{Bmatrix}\frac{\partial{\Psi_{I}^{e\; 1}\left( {\xi_{gi},\eta_{{gi}\;}} \right)}}{\partial\xi} \\\frac{\partial{\Psi_{I}^{e\; 1}\left( {\xi_{gi},\eta_{gi}} \right)}}{\partial\eta}\end{Bmatrix}}},{i = {2\mspace{14mu}{and}\mspace{14mu} 3}}} & (68) \\{{\begin{Bmatrix}\left( b_{Ix} \right)_{i}^{e\; 2} \\\left( b_{Ix} \right)_{i}^{e\; 2}\end{Bmatrix} = {J_{i}^{e\; 2^{- 1}}\begin{Bmatrix}\frac{\partial{\Psi_{I}^{e\; 2}\left( {\xi_{gi},\eta_{gi}} \right)}}{\partial\xi} \\\frac{\partial{\Psi_{I}^{e\; 2}\left( {\xi_{gi},\eta_{gi}} \right)}}{\partial\eta}\end{Bmatrix}}},{i = {1\mspace{14mu}{and}\mspace{14mu} 4}}} & (69)\end{matrix}$

It is easy to verify that for multiple elements nl as shown in FIG. 4,the smoothed gradient matrix of inter node I satisfies the integrationconstraint, i.e.

$\begin{matrix}{{\int_{\Omega}{{\overset{\sim}{b}}_{I}{\mathbb{d}\Omega}}} = {{\sum\limits_{n = 1}^{nl}{\int_{- 1}^{1}{\int_{- 1}^{1}{{\overset{\sim}{b}}_{I}^{n}{\det\left( J_{\xi}^{n} \right)}{\mathbb{d}\xi}{\mathbb{d}\eta}}}}} = {{\sum\limits_{n = 1}^{nl}{\int_{- 1}^{1}{\int_{- 1}^{1}{b_{I}^{n}{\det\left( J_{\xi}^{n} \right)}{\mathbb{d}\xi}\;{\mathbb{d}\eta}}}}} = 0}}} & (70)\end{matrix}$

Therefore the area-weighted integration preserves the linear exactnessand the resultant ME-FEM formulation is valid in the compressible case,for example, passing the patch test (one of the fundamental requirementsin finite element analysis).

Modified Hu-Washizu Formulation

To introduce the strain smoothing formulation into the Galerkinapproximation, a weak form of the Hu-Washizu variational principle isconsidered and given byδU _(HW)(u,{tilde over (ε)},{tilde over (σ)})=∫_(Ω)δ{tilde over (ε)}^(T)C{tilde over (ε)}dΩ+∫ _(Ω)δ{tilde over (σ)}(∇u−{tilde over (ε)})dΩ−δW_(ext)  (71)where u is the displacement field, {tilde over (ε)} is the assume strainfield and iii is the stress field evaluated from the constitutiveequation in Eq. (3) using the assumed strain field. The term δW_(ext)designates the external virtual work by the body force f and traction ton the natural boundary Γ_(N).

The assumed strain method is adopted to eliminate the second term on theRHS of Eq. (71) by assuming that stress field {tilde over (σ)} isorthogonal to the difference between the displacement gradient ∇u andsmoothed strain filed {tilde over (ε)}. Using Eq. (60), one can findthat

$\begin{matrix}{{\int_{\Omega}{\delta{\overset{\sim}{\sigma}\left( {{\nabla u} - \overset{\sim}{ɛ}} \right)}{\mathbb{d}\Omega}}} = {{\sum\limits_{m = 1}^{n\; m}{\delta\;\overset{\sim}{\sigma}{\int_{\Omega_{m}}{\left( {{\nabla u} - \overset{\sim}{ɛ}} \right){\mathbb{d}\Omega}}}}} = 0}} & (72)\end{matrix}$which exactly satisfies the orthogonal condition. The index nm in Eq.(72) denotes the total number of smoothed domain Ω_(m) in Ω(Ω=∪_(m)Ω_(m)).After eliminating the stress components from Eq. (71) using Eq. (72),the following weak form of the modified Hu-Washizu functional isobtained depending only on strain and displacement fields:δU _(HW mod)(u,{tilde over (ε)})=∫ _(Ω)δ{tilde over (ε)}^(T) C{tildeover (ε)}dΩ−δW _(ext)  (73)The displacements in each element are approximated by the ME-FEMinterpolation and are expressed by

$\begin{matrix}{u^{h} = {\sum\limits_{I = 1}^{5}{\Psi_{I}u_{I}}}} & (74)\end{matrix}$The strains are approximated and given by

$\begin{matrix}{{\overset{\sim}{ɛ}}^{h} = {\sum\limits_{I = 1}^{NP}{{\overset{\sim}{B}}_{I}u_{I}}}} & (75)\end{matrix}$where {tilde over (B)} is the smoothed gradient matrix defined in Eq.(65) and 5≦NP≦8 is the number of nodes involved in each smoothed domainΩ_(m). By introducing the displacement and strain approximations intoEq. (73), the following discrete governing equation is obtained:{tilde over (K)}d=f ^(ext)  (76){tilde over (K)} _(IJ)=∫_(Ω) {tilde over (B)} _(I) ^(T) C{tilde over(B)} _(J) dΩ  (77)f ^(ext)=∫_(Ω)Ψ^(T) fdΩ+∫ _(Γ) _(N) Ψ^(T) tdΓ  (78)

Since ME-FEM element interpolation preserves Kronecker-delta property atthe element boundary, essential boundary conditions can be applied in anormal manner (e.g., finite element analysis).

When the eigenvalue analysis for a single element is reanalyzed usingME-FEM method including strain smoothing, almost identical locking-freemodes shown in FIG. 5B are obtained.

Similarly, for multiple elements the results are presented in FIG. 7.Comparing to the EFG results in FIG. 1B, eight eigenvalues of ME-FEMelements with strain smoothing are bounded in the near-incompressiblelimit and the corresponding deformation modes are non-locking

Symmetry Properties of Shape Function Derivatives

The symmetry properties of the shape function derivatives in thefive-noded ME-FEM element are expressed byΨ_(1,ξ)(ε_(g1),η_(g1))=−Ψ_(2,ξ)(ξ_(g2),η_(g2))=−Ψ_(3,ξ)(ξ_(g3),η_(g3))=Ψ_(4,ξ)(ξ_(g4),η_(g4))Ψ_(1,η)(ε_(g1),η_(g1))=Ψ_(2,η)(ξ_(g2),η_(g2))=−Ψ_(3,η)(ξ_(g3),η_(g3))=−Ψ_(4,η)(ξ_(g4),η_(g4))Ψ_(2,ξ)(ε_(g1),η_(g1))=−Ψ_(1,ξ)(ξ_(g2),η_(g2))=−Ψ_(4,ξ)(ξ_(g3),η_(g3))=Ψ_(3,ξ)(ξ_(g4),η_(g4))Ψ_(4,η)(ε_(g1),η_(g1))=Ψ_(3,η)(ξ_(g2),η_(g2))=−Ψ_(2,η)(ξ_(g3),η_(g3))=−Ψ_(1,η)(ξ_(g4),η_(g4))Ψ_(2,η)(ε_(g1),η_(g1))=Ψ_(1,η)(ξ_(g2),η_(g2))=−Ψ_(4,η)(ξ_(g3),η_(g3))=−Ψ_(3,η)(ξ_(g4),η_(g4))Ψ_(4,ξ)(ε_(g1),η_(g1))=−Ψ_(3,ξ)(ξ_(g2),η_(g2))=−Ψ_(2,ξ)(ξ_(g3),η_(g3))=Ψ_(1,ξ)(ξ_(g4),η_(g4))Ψ_(3,ξ)(ε_(g1),η_(g1))=−Ψ_(4,ξ)(ξ_(g2),η_(g2))=−Ψ_(1,ξ)(ξ_(g3),η_(g3))=Ψ_(2,ξ)(ξ_(g4),η_(g4))Ψ_(3,η)(ε_(g1),η_(g1))=Ψ_(4,η)(ξ_(g2),η_(g2))=−Ψ_(1,η)(ξ_(g3),η_(g3))=−Ψ_(2,η)(ξ_(g4),η_(g4))Ψ_(5,ξ)(ε_(g1),η_(g1))=−Ψ_(5,ξ)(ξ_(g2),η_(g2))=−Ψ_(5,ξ)(ξ_(g3),η_(g3))=Ψ_(5,ξ)(ξ_(g4),η_(g4))Ψ_(5,η)(ε_(g1),η_(g1))=Ψ_(5,η)(ξ_(g2),η_(g2))=−Ψ_(5,η)(ξ_(g3),η_(g3))=−Ψ_(5,η)(ξ_(g4),η_(g4))Ψ_(1,ξ)(ε_(g1),η_(g1))=Ψ_(1,η)(ξ_(g1),η_(g1))Ψ_(2,ξ)(ε_(g1),η_(g1))=Ψ_(4,η)(ξ_(g1),η_(g1))Ψ_(5,ξ)(ε_(g1),η_(g1))=Ψ_(5,η)(ξ_(g1),η_(g1))where ξ_(gi)=(ξ_(gi),η_(gi)) i=1, 4 are Gauss points 610.

Referring now to FIG. 8, it is shown a flowchart illustrating anexemplary process 800 for numerical simulation of structural behaviorsin compressible or near-incompressible region using a meshfree-enrichedfinite element method in accordance with one embodiment of the presentinvention. Process 800 can be implemented in software or an applicationmodule executed in a computer system.

Process 800 starts, at step 802, receiving a FEM model of an engineeringproduct (e.g., car, airplane, or components). The FEM model containslow-order FEM elements, for example, bilinear quadrilateral elements ortriangular elements. These low-order elements contain only corner nodeswithout any mid-edge or middle node.

Next, at step 804, a meshfree enriched FEM (ME-FEM) model is created byadding at least one meshfree enriched (ME) node to each of a subset ofthe FEM elements. The ME node is configured to use generalized meshfreeapproximation to eliminate volumetric locking described in previousparagraphs. The subset can include the entire or just a portion of theFEM model. In other words, a mixed mode of ME-FEM and FEM elements maybe used for the simulation.

Process 800 then performs a strain smoothing procedure usingarea-weighted integration technique to avoid potential pressureoscillation problem (a numerical problem) at step 806. Applying suchtechnique is to ensure the simulation results are more realistic withoutusing postprocessing averaging technique.

Finally, at step 808, after the ME-FEM model has been created,numerically simulated structural behaviors of the engineering productare obtained. This is achieved by using the ME-FEM model with anapplication module installed in a computer system. The applicationmodule is configured for performing a meshfree enriched finite elementanalysis that uses a displacement-based first-order convex meshfreeapproximation and an element-wise interpolation, The interpolationpossesses Kronecker-delta property at boundary of each ME-FEM element toallow essential boundary conditions to be applied. Thedisplacement-based first-order convex meshfree approximation isconfigured to eliminate volumetric locking. The numerically simulatedstructural behaviors can be used for assisting user (e.g., scientist,engineer, etc.) to improve the engineering product.

According to one aspect, the present invention is directed towards oneor more computer systems capable of carrying out the functionalitydescribed herein. An example of a computer system 900 is shown in FIG.9. The computer system 900 includes one or more processors, such asprocessor 904. The processor 904 is connected to a computer systeminternal communication bus 902. Various software embodiments aredescribed in terms of this exemplary computer system. After reading thisdescription, it will become apparent to a person skilled in the relevantart(s) how to implement the invention using other computer systemsand/or computer architectures.

Computer system 900 also includes a main memory 908, preferably randomaccess memory (RAM), and may also include a secondary memory 910. Thesecondary memory 910 may include, for example, one or more hard diskdrives 912 and/or one or more removable storage drives 914, representinga floppy disk drive, a magnetic tape drive, an optical disk drive, etc.The removable storage drive 914 reads from and/or writes to a removablestorage unit 918 in a well-known manner. Removable storage unit 918,represents a floppy disk, magnetic tape, optical disk, etc. which isread by and written to by removable storage drive 914. As will beappreciated, the removable storage unit 918 includes a computer readablestorage medium having stored therein computer software and/or data.

In alternative embodiments, secondary memory 910 may include othersimilar means for allowing computer programs or other instructions to beloaded into computer system 900. Such means may include, for example, aremovable storage unit 922 and an interface 920. Examples of such mayinclude a program cartridge and cartridge interface (such as that foundin video game devices), a removable memory chip (such as an ErasableProgrammable Read-Only Memory (EPROM), Universal Serial Bus (USB) flashmemory, or PROM) and associated socket, and other removable storageunits 922 and interfaces 920 which allow software and data to betransferred from the removable storage unit 922 to computer system 900.In general, Computer system 900 is controlled and coordinated byoperating system (OS) software, which performs tasks such as processscheduling, memory management, networking and I/O services.

There may also be a communications interface 924 connecting to the bus902. Communications interface 924 allows software and data to betransferred between computer system 900 and external devices. Examplesof communications interface 924 may include a modem, a network interface(such as an Ethernet card), a communications port, a Personal ComputerMemory Card International Association (PCMCIA) slot and card, etc. Thecomputer 900 communicates with other computing devices over a datanetwork based on a special set of rules (i.e., a protocol). One of thecommon protocols is TCP/IP (Transmission Control Protocol/InternetProtocol) commonly used in the Internet. In general, the communicationinterface 924 manages the assembling of a data file into smaller packetsthat are transmitted over the data network or reassembles receivedpackets into the original data file. In addition, the communicationinterface 924 handles the address part of each packet so that it gets tothe right destination or intercepts packets destined for the computer900. In this document, the terms “computer program medium” and “computerreadable medium” are used to generally refer to media such as removablestorage drive 914, and/or a hard disk installed in hard disk drive 912.These computer program products are means for providing software tocomputer system 900. The invention is directed to such computer programproducts.

The computer system 900 may also include an input/output (I/O) interface930, which provides the computer system 900 to access monitor, keyboard,mouse, printer, scanner, plotter, and alike.

Computer programs (also called computer control logic) are stored asapplication modules 906 in main memory 908 and/or secondary memory 910.Computer programs may also be received via communications interface 924.Such computer programs, when executed, enable the computer system 900 toperform the features of the present invention as discussed herein. Inparticular, the computer programs, when executed, enable the processor904 to perform features of the present invention. Accordingly, suchcomputer programs represent controllers of the computer system 900.

In an embodiment where the invention is implemented using software, thesoftware may be stored in a computer program product and loaded intocomputer system 900 using removable storage drive 914, hard drive 912,or communications interface 924. The application module 906, whenexecuted by the processor 904, causes the processor 904 to perform thefunctions of the invention as described herein.

The main memory 908 may be loaded with one or more application modules906 that can be executed by one or more processors 904 with or without auser input through the I/O interface 930 to achieve desired tasks. Inoperation, when at least one processor 904 executes one of theapplication modules 906, the results are computed and stored in thesecondary memory 910 (i.e., hard disk drive 912). The status of themeshfree enriched finite element analysis is reported to the user viathe I/O interface 930 either in a text or in a graphical representationto a monitor coupling to the computer system 900.

Although the present invention has been described with reference tospecific embodiments thereof, these embodiments are merely illustrative,and not restrictive of, the present invention. Various modifications orchanges to the specifically disclosed exemplary embodiments will besuggested to persons skilled in the art. For example, whereas afive-noded quadrilateral element (i.e., an element having four (4)corner nodes plus one (1) meshfree enriched node) has been shown anddescribed to explain how ME-FEM works mathematically. Other types ofME-FEM element can be used for achieving substantially similarobjectives of the present invention, for example, six-nodedquadrilateral element (i.e., element having four (4) corner nodes plustwo (2) ME nodes), or a four-noded triangular element. Additionally,whereas two-dimensional element (quadrilateral element) has been shownand described to demonstrate how the present invention works, thepresent invention can be implemented on a three dimensional element(e.g., tetrahedral element shown in FIG. 2D). ME-FEM for triangular andtetrahedral elements are substantially similar to that of quadrilateralelements. In summary, the scope of the invention should not berestricted to the specific exemplary embodiments disclosed herein, andall modifications that are readily suggested to those of ordinary skillin the art should be included within the spirit and purview of thisapplication and scope of the appended claims.

We claim:
 1. A computer-implemented method of simulating structuralbehaviors of an engineering product, the method comprising: receiving,by an application module installed in a computer system, a finiteelement method (FEM) model representing an engineering product, said FEMmodel comprises a plurality of FEM elements; creating, by saidapplication module, a meshfree enriched FEM (ME-FEM) model by adding atleast one meshfree enriched (ME) node into each of a subset of the FEMelements, said at least one ME node is configured to eliminatevolumetric locking; and obtaining, by said application module,numerically-simulated structural behaviors of the engineering productusing the ME-FEM model, wherein said application module is configured toperform a meshfree enriched finite element analysis that employs adisplacement-based first-order convex meshfree approximation and anelement-wise interpolation possessing Kronecker-delta property at saideach of said subset of the FEM element's boundary allowing essentialboundary conditions to be applied.
 2. The computer-implemented method ofclaim 1, wherein the FEM elements comprise low-order finite elementsbeing defined by corner nodes and edges, said edges are defined byconnecting a pair of adjacent said corner nodes.
 3. Thecomputer-implemented method of claim 1, wherein said displacement-basedfirst-order convex meshfree approximation uses a basis function derivedfrom one of an exponential, a hyperbolic tangent, an inverse tangent andRenyi function.
 4. The computer-implemented method of claim 1, whereinsaid at least one ME node is located within corresponding domain of saideach of said subset of the FEM elements.
 5. The computer-implementedmethod of claim 1, wherein said each of said subset of the FEM elementsis mapped to a transformed space to enable applying one or more domainsof influence for the at least one ME node.
 6. The computer-implementedmethod of claim 5, wherein the transformed space is a parametric spacecorresponding to one of said domains of influence.
 7. Thecomputer-implemented method of claim 1, wherein said structuralbehaviors comprises pressures and strains.
 8. The computer-implementedmethod of claim 7, wherein said ME-FEM further uses an area-weightedintegration incorporated with the element-wise interpolation to smoothensaid pressures and said strains.
 9. The computer-implemented method ofclaim 8, wherein the smoothed strains are added into meshfreeformulation using a modified Hu-Washizu variation principle.
 10. Thecomputer-implemented method of claim 8, wherein said element-wiseinterpolation has a discrete divergence-free property.
 11. Anon-transitory computer-readable storage medium containing instructionsfor controlling a computer system to perform a simulation of structuralbehaviors of an engineering product by a method comprising: receiving,by an application module installed in a computer system, a finiteelement method (FEM) model representing an engineering product, said FEMmodel comprises a plurality of FEM elements; creating, by saidapplication module, a meshfree enriched FEM (ME-FEM) model by adding atleast one meshfree enriched (ME) node into each of a subset of the FEMelements, said at least one ME node is configured to eliminatevolumetric locking; and obtaining, by said application module,numerically-simulated structural behaviors of the engineering productusing the ME-FEM model, wherein said application module is configured toperform a meshfree enriched finite element analysis that employs adisplacement-based first-order convex meshfree approximation and anelement-wise interpolation possessing Kronecker-delta property at saideach of said subset of the FEM element's boundary allowing essentialboundary conditions to be applied.
 12. A system for simulatingstructural behaviors of an engineering product, the system comprising: amain memory for storing computer readable code for an applicationmodule; a plurality of processing units coupled to the main memory, saidprocessing units executing the computer readable code in the main memoryto cause the application module to perform operations by a method of:receiving a finite element method (FEM) model representing anengineering product, said FEM model comprises a plurality of FEMelements; creating a meshfree enriched FEM (ME-FEM) model by adding atleast one meshfree enriched (ME) node into each of a subset of the FEMelements, said at least one ME node is configured to eliminatevolumetric locking; and obtaining numerically-simulated structuralbehaviors of the engineering product using the ME-FEM model, whereinsaid application module is configured to perform a meshfree enrichedfinite element analysis that employs a displacement-based first-orderconvex meshfree approximation and an element-wise interpolationpossessing Kronecker-delta property at said each of said subset of theFEM element's boundary allowing essential boundary conditions to beapplied.